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ABSTRACT 


The  shock  layer  about  a  blunt  body  in  a  hypersonic  flow 
field  is  studied  considering  both  viscous  and  chemical  nonequi¬ 
librium  effects.  An  order  of  magnitude  study  is  conducted  to 
predict  flight  regimes  for  chemical  equilibrium  and  for  frozen 
flow.  Conservation  equations  are  formulated  for  the  viscous, 
reacting  flow  in  the  shock  layer  about  a  sphere  using  a  seven 
species,  six  reaction  chemical  model.  A  numerical  method  for 
solution  of  these  equations  is  developed  and  the  results  of 
solution  attempts  discussed. 
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I.  INTRODUCTION 


Conversion  of  the  kinetic  energy  of  a  hypersonic  vehicle  to 
thermal  energy  through  shock  waves  and  viscous  effects  may  result  in 
dissociation  and  ionization  of  the  gas  surrounding  the  vehicle.  A 
knowledge  of  the  properties  and  extent  of  this  ionized  layer  is  important 
in  determining  energy  transfer  to  the  vehicle  and  methods  of  communi¬ 
cation  with  it.  In  order  to  determine  the  structure  of  the  reacting  layer 
one  must  take  intb  account  the  body  shape,  shock  shape,  nonequilibrium 
dissociation  and  recombination,  rotational  and  vibrational  relaxation, 
diffusion  and  viscous  effects.  Most  attempts  at  an  analytic  formulation 
of  the  problem  have  assumed  that  viscosity  effects  are  limited  to  a 
relatively  thin  boundary  layer  within  the  shock  layer.  For  bodies  at  high 
altitudes  this  model  is  inappropriate  as  viscous  effects  are  significant 
throughout  the  entire  shock  layer.  The  purpose  of  this  report  is  to 
determine  the  flight  regimes,  as  specified  by  Mach  number  and  altitude, 
in  which  viscosity  and  nonequilibrium  reaction  effects  are  important  in 
the  shock  layer  about  a  blunt  body,  and  to  attempt  an  analytic  or  numerical 
solution  to  the  flow.  In  order  to  accomplish  this  object  it  is  necessary  to 
determine,  by  a  literature  search,  what  reactions  would  be  important  and 
what  the  reaction  rate  constants  for  these  reactions  should  be. 


A  literature  survey  of  current  work  in  flow  about  blunt  bodies  of 
hypersonic  speeds  revealed  that  most  authors  have  considered  either  an 
inviscid,  nonequilibrium  flow  or  a  viscous  flow  in  chemical  equilibrium. 
The  analytical  models  used  are  for  the  most  part  based  on  a  boundary- 
layer  flow,  with  curvature  and  normal  pressure  gradient  effects  neglected. 
Chung  (Ref.  1)  has  given  solutions  for  flow  in  the  viscous  shock  layer  with 
nonequilibrium  chemical  reactions  for  an  "air  atom"  and  "air  molecule" 
model.  His  analysis  is  based  on  a  thin  shock  layer  with  boundary -layer 
type  assumptions.  A  solution  for  nonreacting  flow  using  the  complete 
Navier-Stokes  equations  has  been  obtained  by  Ho  and  Probstein  (Ref.  2) 
by  assuming  forms  for  the  dependent  variables  in  the  neighborhood  of 
the  stagnation  point  such  that  the  Navier-Stokes  equations  reduce  to  a 
set  of  ordinary  differential  equations. 

A  more  recent  review  by  Cheng  (Ref.  3)  of  hypersonic  flow  research 
indicates  that  no  solutions  are  yet  available  in  which  viscous  and  chemical 
nonequilibrium  effects  are  accounted  for  throughout  the  shock-layer 
region.  Additional  references  to  the  literature  will  not  be  presented  here 
since  the  comprehensive  survey  by  Cheng  (ibid)  is  currently  available. 
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II.  FORMULATION  OF  THE  PROBLEM 


A.  A  SIMPLE  ORDER  OF  MAGNITUDE  STUDY 

Before  attempting  to  solve  the  complete  problem  of  the  chemically 
nonequilibrium  viscous  shock  layer  of  the  reentry  body,  a  simple  order  of 
magnitude  study  based  on  the  characteristic  dissociation  and  convective 
times  has  been  performed  (Ref. 4).  From  this  preliminary  analysis  it  is 
possible  to  construct  a  flight  spectrum  to  indicate  the  regions  in  the 
altitude-Mach  number  plane  where  the  chemical  state  of  the  gas  at  the 
stagnation  region  is  expected  to  be  in  a  near  frozen  or  in  a  near  equilibrium 
state.  This  approximate  analysis  is  somewhat  similar  to  those  of  Pearce 
et  al  (Ref.  5)  and  Grier  and  Sands  (Ref.  6).  These  studies  deal  with  the 
recombination  of  species  in  the  boundary  layer  and  are  valid  at  lower 
altitudes  where  the  shock  layer  is  inviscid.  In  the  present  analysis  we 
consider  an  altitude  range  where  the  shock  layer  is  not  inviscid.  The 
entire  shock  layer  is  viscous  and  the  chemical  state  of  the  gas  is  governed 
by  the  rate  of  dissociation.  The  following  assumptions  are  made  in  order 
to  estimate  the  characteristic  times: 

(a)  The  shock  wave  is  a  thin  discontinuity  and  no  reaction  takes 
place  within  the  shock  thickness. 

(b)  The  temperature  of  the  body,  Tg,  is  constant  (1500°K)  and  the 
body  is  noncatalytic. 
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(c)  The  average  thermodynamic  properties  in  the  shock  layer 
determine  the  chemical  state  of  the  gas. 

(d)  In  the  shock  layer  O2  dissociates  first. 

(e)  The  shock  layer  is  essentially  frozen  if  less  than  10% 
dissociation  takes  place. 

Because  of  the  thin  discontinuity,  assumption  (a)  is  self-explanatory. 
A  typical  temperature  of  a  cooled  body  (may  be  ablative)  is  given  by  (b). 
Assumption  (c)  would  be  valid  if  proper  averages  are  taken.  Unfortunately, 
the  correct  average  values  may  not  be  known  accurately  a  priori,  but  a 
reasonable  estimate  for  the  viscous  shock  layer  can  be  made  as  a  zeroth 
approximation.  The  velocity  distribution  along  the  stagnation  streamline 
is  roughly  linear  (us  behind  the  shock  and  zero  at  the  stagnation  region). 

Therefore  u  —  — .  If,  however,  we  use  T  ~  — — 5_,  we  find  that  the 
2  2 

resulting  T  is  lower  than  that  of  the  actual  inviscid  value  (the  subscript  s 
denotes  the  state  just  behind  the  shock).  Because  of  heat  conduction,  the 
average  temperature  of  the  viscous  shock  layer,  in  general,  will  be  lower 
than  the  inviscid  case.  Moreover,  due  to  radiative  heat  transfer,  the 

_  Ts  +  Tb 

temperature  may  decrease  further.  Thus  T  ~  - - -  will  be 

reasonably  close  to  the  actual  average  value  in  the  viscous  shock  layer 
with  radiative  heat  transfer. 

The  upper  limiting  case  is  also  considered  where  T  ~  Ts.  This 
is  the  worst  possible  condition  where  the  shock  layer  is  subjected  to  the 
maximum  possible  temperature. 
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For  a  given  M  ;  — 2.,  ,  — — ,  etc.  can  be  found  easily. 

Too  Poo  Poo  Uoo 

From  these  values  u,  T  and  'p  corresponding  to  T  can  be  evaluated. 


During  its  travel  through  the  shock  layer,  the  fluid  will  tend  to  dissociate 


at  the  above  average  conditions. 

Assumption  (d)  is  quite  accurate  for  thermodynamic  calculations 
but  will  not  convey  the  correct  physical  picture  when  one  is  interested  in 
the  formation  of  NO  and  its  subsequent  ionization,  NO  ;=*  NO+  +  e”  . 
However,  for  the  assumed  model,  if  O2  dissociation  is  frozen  to  start 
with,  the  above  ionization  will  not  take  place.  Since  in  the  present 
simplified  analysis  we  want  to  know  whether  the  gas  is  in  a  near  frozen 
chemical  state  or  in  a  state  of  nonequilibrium,  assumption  (d)  is  consistent 
with  the  other  assumptions. 

Instead  of  using  10%  O^  dissociation  as  the  criterion  of  a  near 
frozen  state, ^assumption  (e),j  we  can  also  assume  1%  dissociation  as  the 
criterion  in  order  to  be  on  the  conservative  side.  That  is  why  Fig.  2 


shows  both  1%  and  10%  O2  dissociation  lines. 

Using  the  relationship  ~  — —  (A  =  shock  standoff  distance 

R-B  Ps 

and  Rg  =  body  radius)  we  can  estimate  the  convective  time  tc  as  . 

u 

Thus  for  a  one-step  reaction  £  v\  M;  M,,  the  ratio  of 

j  J  J  kb  j  J  J 

td 

the  characteristic  dissociation  and  the  convective  times, - can  be 

*c 

given  by  the  following  expression 
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a  =  speed  of  sound,  ft/sec 
M  =  Mach  number  -  ~  ‘ 

For  a  very  large  p  we  get  a  near  frozen  state  and  for  a  very  small  p  we 
get  a  near  equilibrium  state  of  the  shock  layer. 

Eschenroeder,  et  al  (Ref.  7),  give  the  following  rate  constant  for 


02  dissociation. 


o2  +  N2 


20  +  N2.  kf  =  11.9  x  10' 


.-3/2 


59,380> 


mole-sec 


It  is  assumed  here  that  this  expression  can  be  extrapolated  to  temperatures 
as  high  as  15,000°K  or  higher. 

Figure  1  is  a  plot  of  P vs  M  ^  for  various  altitudes.  Figure  2  is 
the  flight  spectrum  cross  plotted  from  Fig.  1.  These  are  for  Rg  =  1  ft.  In 
order  to  be  on  the  conservative  side  the  highest  available  rate  constant  for 
02  dissociation  is  chosen.  The  expression  used  for  the  shock  standoff 
distance  also  gives  somewhat  high  values.  Thus  the  resulting  tc  will  be 
larger  than  the  actual  value.  The  combination  of  these  conservative 


characteristic  times  gives  higher  values  of  percent  02  dissociation. 

_  Ts  +  Tb 

Based  on  an  average  temperature  T  **  -  the  analysis 

•  2 

shows  that  02  dissociation  will  be  nearly  frozen  at  higher  altitudes 


particularly  in  the  viscous -layer  regime.  The  same  basic  conclusion 
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Fig.  1 


56  MILES 


MoO 


iOti  Parameter  For  Rg  =  1  Foot. 
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TITUDE,  MILES 


Fig.  2  Chemical  State  Of  The  Shock  Layer  At  The  Stagnation 
Region  Of  A  1-Foot  Radius  Noncatalytic  Spherical  Body 
Tg  =  1500° K,  ARDC  Model  Atmosphere  (1959). 
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is  also  obtained  by  considering  T  ~  T8.  This  case  is  the  upper  limit 
where  the  average  gas  temperature  in  the  shock  layer  is  equal  to  Ts. 
Obviously  this  upper  limit  is  physically  impossible.  The  results 
essentially  indicate  that  in  the  viscous-layer  regime  ionization  at  the 
stagnation  region  of  a  hypersonic  1-foot  radius  reentry  body  is  unlikely. 

In  view  of  the  extremely  high  speed  of  the  body  this  conclusion  is  contrary 
to  popular  belief.  However,  recent  shock  tube  experiments  of  Herzberg 
(Ref.  8)  indicated  that  the  impurities  such  as  diaphragm  materials  are  the 
main  sources  of  ionization.  By  carefully  controlling  the  impurities  he  was 
able  to  decrease  the  degree  of  ionization  to  a  level  where  the  instruments 
-tabled  to  detect  any  ionized  species.  The  results  obtained  by  Chung  (Ref.  1) 
also  indicate  that  Fig.  2  gives  a  qualitatively  correct  picture.  Two  points 
from  his  study  are  plotted  in  Fig.  2. 

The  effect  of  the  presence  of  atomic  oxygen  in  the  free  stream  is 
also  considered.  Assuming  that  the  atomic  oxygen  is  50  percent  ol  the 
total  oxygen  in  the  free  stream,  the  same  basic  conclusion  is  reached. 

An  order  of  magnitude  study  such  as  ours  is  inherently  approximate. 
But  it  shows  the  correct  trend.  At  lower  altitudes  the  results  will  be  in 
considerable  error.  This  is  because  chemical  equilibrium  will  prevail  in 
the  inviscid  portion  of  the  shock  layer.  For  this  reason,  the  line  showing 
100%  dissociation  will  not  portray  an  accurate  state  of  the  gas. 
Possibly  dissociation  will  also  be  significant  in  this  region. 
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B.  CONSERVATION  EQUATIONS 


While  the  result  of  the  previous  section  gives  the  qualitative  chemical 
nature  of  the  stagnation  region  of  a  l— foot  radius  body,  it  fails  to  give  any 
information  on  the  structure  of  the  shock  layer.  It  also  fails  to  take  into 
account  the  important  chemical  reactions  other  than  the  dissociation  of  O^. 

The  standard  conservation  equations  for  a  reacting  mixture  are 
written  in  spherical  coordinates  for  an  axisymmetric  spherical  body  shown 
in  Fig.  3.  The  assumptions  which  led  to  the  given  form  of  a  particular 
equation  are  listed  along  with  each  equation  (steady  state  system). 

(i)  Conservation  of  overall  mass: 


r2 


A 

dr 


(p  v?‘  u)  + 


1 


r  sinO 


jL 

de 


(vp  sin6) 


0 


Fig.  3  Coordinate  System. 
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where  p  =  mass  density 

~u  =  r-component  of  the  velocity 
v  =  6-component  of  the  velocity. 

(ii)  Conservation  of  individual  species  (assume  uncoupled 
concentration  diffusion  only): 

jz  X  (p  Cl  r2“l +  T^i  is  (p  c‘v  ,i"e) 

.  J_  J.  (JL.  r2  Jzi)  .  _i -  4-  ( — -  —V  TMiI 

r2  dr  V  Sc  dr  i  r  sin  6  ^0  \  Sc  r  de/  *• 

where  c^  =  mass  fraction  of  species  i 

Sc  =  Schmidt  "number  . ..  . 

jjl  =  viscosity  coefficient 

[Mi]  =  rate  of  change  of  the  molal  concentration 
of  species  i 

Mi  =  molecular  weight  of  species  i. 

The  expression  for  the  volumetric  production  rate  of  species  i,  [Mi]  can 
be  written  only  when  the  reaction  equations  are  specified.  For  [mJ»  one 
essentially  writes  the  law  of  mass  action  for  a  given  set  of  reaction 
equations  which  affect  Mi.  These  equations  are  given  later  in  another 
section. 


M 
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(iii)  Conservation  of  r-momentum  (assume  linear  stress -strain 
relationship,  negligible  second  viscosity  coefficient): 


where  p  =  pressure. 


(iv)  Conservation  of  8-momentum  (assume  linear  stress-strain 
relationship,  negligible  second  viscosity  coefficient): 


♦7  '£(*)♦ 

+  7  is  [-p  *  *(r  -35  +  7)] 


l  du 

7  3e 


}] 


cote}  +  3^{r-^(^)+i-|H} 


(v)  Conservation  of  energy  (assume  diffusive  enthalpy  flux  is 
small  compared  to  the  convective  flux  and  no  radiation  heat 
transfer): 

fr) 

+  TtfaW  "$s(  r  "Hr  *1”8)  +  ZOre.etton+u^  +  fll+4' 
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where  k  =  thermal  conductivity 

cp  =  mean  specific  heat  of  the  mixture  at 
constant  pressure  =  £  cicPi 

^reaction  =  net  heat  released  due  to  reaction 
<J>  =  dissipation  function. 


The  expression  for  ^Qreaction  can  he  written  for  a  given  set  of  reaction 
equations.  Thi6  will  be  given  later  after  the  particular  reaction  equations 
are  specified.  The  dissipation  function  <f>  is  given  by 


$  =  P 


H 


dn\z 

*37 ) 


( l  dv 

*17*30 


v  cot  0 


C.  REACTION  MODEL  AND  RATE  EXPRESSIONS 

Atallah  (Ref.  9)  has  surveyed  the  existing  literature  on  the  chemical 
kinetics  of  air  after  a  hypersonic  shock.  He  gives  the  important  reactions 
and  their  appropriate  rate  expressions  based  on  the  best  known  exper¬ 
imental  data.  There  seems  to  be  a  wide  agreement  in  the  literature  as  to 
which  chemical  species  should  be  abundant  due  to  a  hypersonic  shock  in 
air.  However,  a  great  deal  of  uncertainty  exists  in  the  expressions  for 
the  rate  constants.  In  general,  the  experimental  values  are  applicable 
within  a  temperature  range  of  2000°K  to  6000°K.  The  temperature  after 

the  shock  in  front  of  the  reentry  body  may  exceed  the  upper  limit  by  a  large 
margin.  Thus,  one  iB  forced  to  use  the  rate  constants  in  a  temperature 
range  where  their  validity  is  questionable.  Any  error  in  the  reaction  rate 
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will  introduce  a  considerable  amount  of  error  in  the  solution  of  the  shock 


layer.  Bortner  (Ref.  10)  discusses  this  very  point.  Nawrocki  (Ref.  11) 
among  others  gives  an  excellent  review  of  the  state  of  knowledge  of 
reaction  rates.  Since  reviews  such  as  this  are  already  available  in  the 
published  literature,  no  attempt  will  be  made  here  to  survey  the  existing 
experimental  data. 

The  following  six  chemical  reactions  involving  seven  different 
species  are  assumed  to  play  an  important  part  in  the  behavior  of  the 
shock  layer. 


02  + 

N25=*  20  +  N2 

(1) 

O  + 

N25=*  O  +  2N 

(2) 

n2  + 

Os=*NO+  N 

(3) 

o2  + 

N  «==*  NO  +  O 

(4) 

2N  + 

O  «=*  NO  +  N 

(5) 

N  + 

O  NO+  +  e" 

(6) 

Ionization  reactions 

O  «=*  0+  +  e"  and 

N  5=*N' 

+  e”  have  not  been  in¬ 


cluded  in  the  above  set  of  reaction  equations.  These  two  reactions  take 
place  at  extremely  high  temperatures  and  are  important  in  calculating 
transport  coefficients,  specific  heats,  etc.  But  since  we  are  interested 
in  finding  out  if  we  have  ionization  in  the  shock  layer  our  investigation 
should  be  centered  around  NO  ionization,  NO  «=sNO+  +  e".  This  reaction 
occurs  at  a  relatively  lower  temperature  and  is  not  as  vital  for 

thermodynamic  calculations  as  for  locating  the  position  of  antenna  for 
c  ommunic  ation. 
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Having  established  the  pertinent  reaction  equations  we  are  in  a 
position  to  write  expressions  for  [Mi]  for  each  of  the  seven  species. 
Using  the  notations  of  the  previous  section  and  applying  the  law  of  mass 
action  for  each  of  the  six  equations  we  get  the  following  results: 

^  *  -**1 HN + kb,[°]2H  -  *«[°2In]+ “wMM 

=  2kfl  MM  -  2kb![o]2M  -  kf3  [N2][o]+kb3[NO][N] 

+  kf4  [02][n]  -  kb4  [NO]  [O]  -  kf3[N]2[o]  +  kb5  [NO][N] 

"  kf6  [N][°]  +  kb6  [N°+][e“] 

^  =  ‘“fa  [°][Nj]  +  »b2[0][N]2  -  kj,  [Nj[0]  +  kb3  [NO][N] 

^p-=  2kf2[0][N2]  -  2kb2[0][N]2  +  kf3  [N2][0]  -  kb3  [n<5|[ n] 

-kf4  MM  +  kb4  [NO][o]  -  kf5  [0][N]2  +  kb5  [NO][N] 
-kf6[N][°]+  k&6  [NO+]  [e“] 

=  kf3  [N2][0]  >kb3[NO][N]+  kf4[02][N]-kb4[N0][6] 

+  kf5[N]2  [O]  -  kbs  [no][n] 
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=  kf6  [n][°]  -  kb6[NO+][.-] 

f1 

The  subscripts  f  and  b  refer  to  the  forward  and  the  reverse  reaction  rate 
constants  respectively.  The  number  in  the  subscript  refers  to  the  reaction 
equation.  For  example,  kf  j  refers  to  the  forward  reaction  rate  constant 
for  reaction  3. 

Now  ]£  Qreaction  *n  *ke  ePergy  equation  is  the  net  heat  released 
or  absorbed  due  to  reaction.  It  can  be  evaluated  either  by  summing  up  all 
the  heats  of  reaction  for  individual  reaction  products  or  by  summing  up 
each  molecular  species  multiplied  by  its  heat  of  formation. 

^reaction  •  -  {  1  [°2][N2]  -  kb,  [O]2  [N2]}  Ah, 

-{*f2[°]K]  -  M0]^}  ^2 

-{kf3[N£|[Oj  -  kb3[Nq[N]}Ah3  +{ki4  [0^[Nj-  kb4[NO][0]}Ah4 
+  {kf5[Nf[0]  -  kb5[NO][N]}Ah5  -{k£6[N][0]-kb6[N0+J[e-]}Ah6 

where  Ah's  are  molal  heats  of  reaction. 

The  other  way  of  expressing  the  net  heat  release  is  £  ^reaction  = 

X[Mi]  h*  ;  where  hi  is  the  molal  heat  of  formation  of  species  i.  Because 
of  the  definition  of  the  heat  of  formation  only  the  molecular  species,  O2 
and  N2  need  not  be  considered. 


*  1 6- 


Bortner  (Ref.  12)  gives  the  following  heats  of  reaction  for  the  six 


reaction  equations  we  have: 


°2  +  n2  — -  20  +  N2 
O  +  N2  —  O  +  2N 
N2  +  O  —  NO  +  N 
02  +  N  — — NO  +  O 
2N  +  O  — "NO  +  N 
N  +  O  —  NO+  +  e" 


.  kg -cal 
Ah*  g-mole 

117.7 

226.0 

76.1 

-  32.3 

-  150.0 

64.5 


-Ah  signifies  an  exothermic  reaction. 


The  rate  constants  for  all  the  reactions  except  the  fifth  reaction  are 
taken  from  Ref.  9-  The  rate  constants  for  the  fifth  reaction  are  taken 
from  AVCO  Report  104. 


kfl 


2.36  x  1020  T"3/2  exp( 


59,380  \ 
T  / 


cc 

mole- sec 


kb  =  1.86  x  1017  T”1  - — - 

U1  2 

mole  -sec 


20  —-3/2  /  1 13,340  \  cc 

kf2  =  7.35  x  !02°  T  V2  .*p(-  — J— )  -J..— 


kb?  =  3.9  x 


1019  I'3/2  - SSl 


mole  -sec 
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kf.  =  5.5  x  10“  T  '  exp 


11  t1/2  /  36iM5\  cc 


(.3^45) 


T  /  mole-sec 


kb,  =  2x  1011  T1/2  — ^ - 

3  mole-sec 


*i4  =  8.9  x  10U  T1/2  «p(-^) 


CC 


mole -sec 


kb4  =  1.21  xlO11  T1/2  exp^  - 


750  \  cc 


T  /  mole-sec 


kf  =  2.0  x  10fci  T 


21  T -3/2 


cc‘ 


mole  -sec 


?i  — ^ /  75,000  \ 

kb5  =  7.9  x  1021  T3/2  exp  -  ) 


cc 


mole-sec 


kf  =  6.4  x  10v  T 


i/2  e,p  (.1M55) 


36,650\  cc 


mole -sec 


kb  =  1.78  x  1019  T_1  - — 


mole -Bee 


The  temperature  T  is  in  degrees  Kelvin. 


D.  AUXILIARY  EQUATIONS  AND  BOUNDARY  CONDITIONS 

The  conservation  equations  listed  in  Section  B  do  not  constitute  a 
soluble  set  unless  other  auxiliary  equations  are  known.  In  this  section, 
the  other  necessary  equations  are  listed. 
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(a)  Equation  of  state: 


Perfect  gas  for  individual  species  pi  =  pci  jry  T 

IMi 


Perfect  gas  for  the  ambient  air  Pqq  = 

Mqo 


Real  gas  for  the  shock  layer  p  =  pZ  T  =  p  T 


where  Z  =  compressibility  factor, 


tfe  x  -pU 

Moo 


M=  mean  molecular  weight  V  — i _ 


(b)  Transport  coefficients  and  properties: 

In  the  absence  of  detailed  knowledge  of  the  transport  coefficient 
in  a  nonequilibrium  chemical  system  the  following  assumptions 
are  made  as  a  zeroth  approximation 

[j,  T  k  T 

~c~  =  “ — 5  r —  =  ^ — ;  Pr  =  1;  Sc  =  lj  Cp  =  const. 

^OO  1  00  Koo  1  00 

For  a  gaseous  mixture  in  chemical  equilibrium,  reasonably 
accurate  estimates  of  the  transport  coefficients  are  available 
in  the  literature  {NACA  TN  4150,  for  example).  It  is  quite 
possible  that  the  above  assumptions  of  the  behavior  of  the 
transport  coefficients  do  not  correctly  describe  our  system. 
However,  their  choice  is  dictated  by  their  simpler  form. 

Other  available  relationships  will  neither  provide  more 
accuracy  nor  will  they  help  us  better  understand  the 
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problem.  If  our  solution  shows  a  trend  toward  near 
equilibrium  chemical  state  then  we  plan  to  use  appro¬ 
priate  known  transport  properties. 

(c)  The  state  behind  the  shock: 

T s*  Ps*  Ps»  Mb»  etc.  can  be  found  from  the  shock  tables  for 
a  given  T^,  p^,  IM  QO,  etc.  We  must  remember  in  this  con¬ 
nection  that  no  reaction  takes  place  in  the  shock  thickness 
(thin  discontinuity).  Thus,  we  have  air  at  a  very  high 
temperature  at  r  =  RB. 

(d)  Shock  radius,  Rg: 

We  shall  assume  ~  0.1.  Therefore,  Rs  ~  Rfi 

+  A  ~  1.1  Rg  where  A  =  shock  standoff  distance  at  the 
stagnation  region.  We  shall  also  assume  that  the  shock  is 
spherical  for  small  0.  Choosing  A  arbitrarily  eliminates 
the  free-boundary  nature  of  the  problem  and  consequently 
the  number  of  boundary  conditions  which  may  be  satisfied  by 
the  solution  will  be  reduced  by  one.  If  in  the  solution  of  the 
equations,  it  appears  impossible  to  satisfy  the  specified 
boundary  conditions  at  the  body  using  A  ~  0.1  Rb  this 
assumption  may  be  varied. 
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(e)  Boundary  conditions: 


(i)  Noncatalytic  body: 

(At  r  =  Rg  and  any  0);  u  =  0;  v  =  0;  T  =  Tg; 


dc{ 


=  0. 


(At  r  =  Rs  and  any  6);  u  =  -u8;  v  =  vs;  T  =  Ts;  p  =  p8 


[ri(7)+7lf] 

r  =  R0 


(no  shearing  stress 
at  the  shock) 


h/  =  °- 


atomic 

species 


Also 


dci 

de 


=  o 

r  =  Rg 


(ii)  Catalytic  body: 

(At  r  =  Rg  and  any  0);  u  =  0;  v  =  0;  T  =  Tg 

c^  other  than  cq2  anc^  CN2  are  zero  (complete 

rec  ombination) 

At  r  =  Rs  the  boundary  conditions  are  the  same  as  (i). 


From  the  set  of  equations  shown  in  the  preceding  sections  we 
find  that  we  have  two  independent  variables,  r  and  0.  The 
numbers  of  dependent  variables  and  equations  are  as  follows: 
Dependent  variables  u,  v,  p,  p,  T,  c^  =  12 

5  7 
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Number  of  equations: 


Conservation  of  mass  (overall)  1 

Conservation  of  species  7 

Conservation  of  momentum  2 

Conservation  of  energy  1 

Equation  of  state  _ 1_ 

Total  12 


Thus,  we  have  twelve  equations  and  twelve  unknowns,  and 
the  system  of  equations  is  soluble. 

Three  altitudes,  175,000  feet,  250,000  feet  and  295,000  feet  and 
three  (Moo,  20,  25  and  30  are  to  be  investigated.  These  altitudes  will 
span  the  "viscous -layer  regime"  and  come  down  to  the  fringe  of  the 
"boundary-layer  regime"  as  defined  by  Probstein  (Ref.  13). 

An  analytical  solution  of  these  equations  was  attempted  by  a  method 
similar  to  that  used  by  Probstein  and  Kemp  (Ref.  14)  in  their  solution  for 
nonreacting,  constant  density  viscous  flow  in  the  shock  layer  about  a 
sphere.  Functional  forms  for  the  dependent  variables  were  assumed  in  an 
attempt  to  reduce  the  partial  differential  equations  to  a  set  of  simultaneous 
ordinary  differential  equations  valid  in  the  stagnation  region.  Both  this 
method  and  a  series  solution  method  involving  expansion  in  6  about  the 
stagnation  line  were  unsuccessful. 

Since  no  similarity  parameters  valid  beyond  the  stagnation  point 
are  available,  a  relaxation  scheme  for  the  shock  layer  has  been  developed 
in  the  next  chapter.  The  equations  we  have  are  elliptic  in  nature  and  are 
suitable  for  relaxation  solution. 
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III.  NUMERICAL  PROCEDURE 


The  boundary- value  problem  consists  of  twelve  nonlinear  partial 
differential  equations  of  the  second  order,  with  twelve  dependent  variables 
and  two  independent  ones.  These  may  be  written  in  the  form: 

pk(yi*  Vv  virr*  Vi6*  vi00*  Vir0'  r»  6  >  =  0  (?) 

i  =  1,  2  ...  12 
k  =  1,  2.  .  .12 

with  boundary  values  at  the  shock  (r  =  Rs)  and  the  body  (r  =  Rg). 

In  general,  the  approach  to  this  problem  is  based  on  two  ideas; 
first,  that  the  differential  equations  can  be  replaced  by  finite  difference 
equations  which  are  algebraic  in  the  dependent  variables  at  a  point;  and 
second,  that  systems  of  simultaneous  nonlinear  algebraic  equations  can 
be  solved  by  Newton's  Method.  Thus, a  solution  consists  of  sets  of  the 
dependent  variables  which  satisfy  Eq.  (7)  at  a  finite  number  of  points  (r,  0) 
between  the  shock  and  the  body,  with  spacing  Ar,  A0.  If  a  solution  also 
satisfies  boundary  values  at  the  body,  it  is  assumed  to  be  a  physical 
solution  of  the  differential  equations. 
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A.  DIFFERENCE  EQUATIONS 


All  partial  derivatives  are  replaced  by  their  corresponding  forward 
difference  quotients; 


<J  vi(r,6)  -  vi(r-Ar,  8) 

dr  Vi  =  Ar 


dZ  •fe(vj(r»e))  (vj(r-Ar,9)) 

d?-  Vi  “  Ar 


d 

He 


vi(r,e)  -  vi(r,  e<-Ae) 
Ae 


-^(vi(r*e0  "-^(Vi(r'0-Ae)) 

ae2  Vi  Ae 


a 

Ur 


(vj(r,6))  -  (  v^r.e-Ae)) 


ar  <3e 


Ae 


If  Vj(r-Ar,6)  and  vj(r,  6-A6)  are  known,  the  problem  is  reduced 
to  twelve  nonlinear  simultaneous  algebraic  equations  in  twelve  unknowns 
VJ  at  any  point  (r,  6). 


Fk(Vi) 


=  0 

r.e 


(8) 
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B.  NEWTON'S  METHOD 


Assume  that  for  fixed  r,  0  a  set  has  somehow  been  estimated 

which  is  "sufficiently  close"  to  the  solution  set  |v^.  Expand  each  of  the 
Ffc  in  Taylor's  Series  about  the  point  {v^},  neglecting  all  except  the  first 
order  terms.  We  have  twelve  linear  equations 


Fk(vij  +  Avi) 


=  Fk<vii)  + 


Avj  + 


dFk 

+  3^1  Av'2 


Assume  that^vij  +  Av^  is  a  solution  set,  so  that  Fj^vj  +  AvJ  =  0. 
Then  in  matrix  form 


•Fltvii) 

• 

1 

1 

1 

dFi 

dFx 

Avj 

1 

1 

1 

dvi 

» 

1 

dv12 

1 

1 

1 

1 

“F12(vii> 

1 

1 

«fu 

dF]2 

1 

i 

1 

Av12 

dvl 

dv12 

- 

- 

— 

_  _ 

or 

[-F]  =  [j]  jwj 

Then 

and  the  set 

vi2  =  vi2  +  Avi 
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will  be  an  improved  estimate  if  the  process  is  convergent.  This  procedure 
is  used  repeatedly,  giving 

Vi  =  vj  ,  +  Av<  . 

*n  ln-l  4n-l 

until  all  the  Fk(vin)  are  sufficiently  close  to  zero  simultaneously. 


C.  JACOBIAN  ELEMENTS 


The  144 


££k 

dvi 


of  the  J-matrix  may  be  determined  by  differentiating 


each  Ffc  twelve  times  and  writing  144  equations;  however  it  is  more 
efficient  to  find  these  numerically.  An  arbitrary  increment  Svj  is  chosen 
for  each  vj  and  the  derivatives  found  from  the  difference  quotients: 


a  .  *k<vl'  ••••  vi~l»vi  +  &vi*  vi+l  ••••  vi2)-Fk(vi...vi. . 

w  = - £ - - - 


D.  BOUNDARY  VALUES 


For  any  numerical  solution  of  a  differential  equation  it  is  necessary 
to  replace  a  boundary-value  problem  with  an  equivalent  initial-value 
problem.  Some  of  the  initial  valueB  will  be  chosen  arbitrarily.  Each 
choice  will  give  a  solution  which  may  or  may  not  match  the  known  conditions 
on  the  final  boundary.  The  arbitrary  initial  values  are  then  varied  until  a 
solution  is  found  which  satisfies  all  boundary  conditions. 


•v12) 


-26- 


In  this  problem  the  dependent  variables  v$  are  known  at  the  shock 


from  free-stream  conditions  and  the  oblique  shock  relations.  Thus  all 

d  j  a2 

r=R- 


are  also  known.  To  get  approximate  values 


for  Vi  the  inviscid  equations  are  solved  at  the  shock.  These  values 
can  then  be  substituted  into  the  viscous  equations  (7)  to  find  values  of 

d2 

t  Vi  so  that  Eqs.  (8)  are  actually  satisfied  at  the  shock, 
dr 

Values  on  the  stagnation  line  (0  =  0,  Rg  <  r  <  Rs)  must  also  be 
known.  Originally  these  were  estimated  from  boundary  values  and  an 
attempt  to  improve  the  estimate  was  made.  It  is  now  felt  that  the 
stagnation  line  must  be  solved  as  a  separate  problem  (see  Section  IV 
Results). 


E.  EXTRAPOLATION 

To  start  the  Newton's  Method  iteration  at  an  unknown  point  (r,  6) 
an  initial  estimate  of  vi  and  all  derivatives  is  required.  This  is  done  by 
assuming  constant  second  derivatives  and  using  second  order  extrapolation 
from  known  points  (r  -  Ar,  0)  (r,  0  -  A0) 
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Vi(r,e)  =  Vi(r-Ar.e)  +  v*  (r-Ar,  0)  Ar  +  2  vi(r-Ar,0) 


and  similarly  for  the  0 -derivatives. 


F.  SEQUENCE  OF  OPERATIONS 


For  an  individual  point  (r,  0); 

1.  Extrapolate  initial  estimates  from  points  (r-Ar,  0) 
and  (r,  0-A0). 

2.  Evaluate  F^  (vj)  and  test  for  zero.  If  a  solution  has 
been  found  proceed  to  next  point. 

3.  If  the  Ffc  are  not  zero,  set  up  and  solve  the  Jacobian 
matrix  for  the  Av^.  Modify  the  v^  and  go  to  step  2. 

This  procedure  is  started  from  the  shock  and  carried  along  a  radial  line 
0  =  const,  until  the  body  is  reached,  and  repeated  for  increasing  0  until 
the  supersonic  region  is  reached.  If  the  solution  does  not  match  the 

/  dv{  \ 

boundary  values  along  the  body,  the  slopes  at  the  shock^-^ — J  are  adjusted 
and  another  solution  is  obtained. 


G.  SUMMARY  OF  ASSUMPTIONS  MADE 

1.  That  a  unique  solution  of  the  differential  equations  exists,  and 
that  solutions  of  the  corresponding  difference  equations  will  converge  to 
this  as  net  size  (Ar,  A0)  decreases. 
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2.  That  the  functions  Fjt(vi)  are  well-behaved  enough  to  allow 
linear  representation  by  Taylor  Series  in  some  neighborhood  of  the  |v^. 

3.  That  initial  estimates  of  the  can  be  made  which  will  fall 
within  this  region  of  convergence  so  that  Newton's  Method  can  be  used. 

4.  That  unknown  boundary  values  can  be  estimated  with  enough 
accuracy  to  permit  assumption  3  to  be  made. 


IV.  RESULTS 

The  complete  flow-field  problem  waB  programmed  for  the 
WDPC  7090  Computer.  The  boundary  line  0  =  -  A0,  Rs  >  r  >  Rg  was 
assumed  to  be  known  and  the  solution  was  attempted  along  the  stagnation 
line.  Although  many  different  sets  of  initial  values  were  tried,  none  was 
accurate  enough  to  fall  within  the  region  of  convergence  of  Newton's 
Method.  In  addition,  since  the  Navier-Stokes  equations  in  spherical  coor¬ 
dinates  become  indeterminant  on  the  stagnation  line  and  a  limiting  process 
must  be  resorted  to,  it  was  decided  to  attempt  to  solve  the  stagnation  line 
separately. 

This  has  two  advantages;  first,  a  stagnation  line  solution  would 
serve  as  a  boundary  for  the  flow-field  problem;  second,  it  provides  a 
simpler  way  to  study  the  operation  of  the  numerical  techniques  used.  It 
was  decided  to  further  simplify  the  problem  by  approaching  it  in  parts; 
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that  is,  to  solve  the  frozen  flow  problem  (three  equations,  three  unknowns) 
first  and  then  use  this  solution  as  a  starting  point  for  the  reacting  flow 
stagnation  line  problem. 

A.  EFFECT  OF  ERROR  ACCUMULATION 

A  program  to  solve  for  (u,  p,  T)  in  cgs  units  was  written  for  the 
Honeywell  800  Computer.  This  program  converges  and  gives  solutions  for 
about  one-half  the  distance  from  the  shock  to  the  body.  It  was  felt  that  part 
of  the  difficulty  with  convergence  of  this  program  was  caused  by  the  loss 
of  numerical  significance.  All  computations  are  floating  point,  to  eight 
significant  digits,  but  change  in  the  last  place  of  u,  p,  or  T  affects  the  higher 
order  digits  of  the  functions. 

In  an  effort  to  eliminate  this  difficulty  the  frozen  flow  problem  was 
reprogrammed  in  dimensionless  form.*  This  improved  the  results  in 
degree  but  not  in  kind.  Although  the  Newton's  Method  solution  at  a  point 
is  self-correcting,  errors  accumulate  from  point  to  point,  as  with  any 
marching  procedure.  It  seems  that  the  increasing  error  degrades  the 
quality  of  the  initial  estimate  at  a  point  until  it  finally  falls  outside  the 
region  of  convergence  of  the  Newton's  Method.  Values  of  the  functions 
of  the  initial  estimates  are  plotted  versus  radius  for  a  typical  case  in 


The  nondimensional  equations  for  frozen  flow  on  the  stagnation 
line  are  given  in  the  Appendix. 
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Fig.  4.  Thus,  this  technique  can  always  be  expected  to  fail  to  converge 
at  some  point.  The  usefulness  of  the  technique  depends  on  whether  or  not 
the  body  is  reached  before  convergence  fails. 

To  attempt  to  make  the  solution  reach  the  assumed  body  position, 
many  sets  of  initial  values  of  slopes  of  u,  p,T  at  the  shock  were  tried. 
Sample  plots  are  shown  in  Figs.  5,  6  and  7.  These  figures  also  show  that 
the  effect  of  changing  the  point  spacing  is  small,  indicating  that  Ar  has 
been  chosen  small  enough.  At  points  where  convergence  was  obtained, 
the  procedure  generally  decreased  the  values  of  the  functions  by  an  order 
of  magnitude  with  each  iteration.  Convergence  at  a  typical  point  is  shown 
in  Figs.  8,  9  and  10.  All  figures  shown  are  for  (Mqq  =  20),  altitude  of 
295,000  feet,  and  body  radius  of  1  foot. 

B.  EFFECT  OF  SHOCK  STANDOFF  DISTANCE 

It  was  felt  at  this  point  that  the  main  difficulty  in  reaching  the  body 
where  u  =  0  was  probably  due  to  the  choice  of  a  fixed  value  of  the  shock 
standoff  distance.  The  actual  shock  standoff  distance  is  a  physical  param¬ 
eter  and  it  should  be  obtained  as  a  part  of  the  solution.  For  a  given  body 
traveling  at  a  specified  speed,  the  shock  standoff  distance  is  a  unique 
distance.  However,  it  is  not  easy  to  determine  the  exact  value  of  this 

*Z - 

Free-stream  Mach  number  of  20. 
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distance  for  a  given  body  and  Mach  number  combination.  The  available 

order  of  magnitude  value  of  -=p-  is  not  accurate  enough  for  a  problem  of 

KB 

this  kind. 

In  our  problem,  one  way  of  determining  this  distance  as  a  part  of 
the  solution  will  be  to  start  the  iteration  from  the  shock  and  proceed 
toward  the  body  through  the  shock  layer.  When  the  solution  gives  us 
u  =  0  and  matches  the  rest  of  the  physical  boundary  conditions  we  have 
actually  reached  the  body.  Thus,  we  get  both  the  structure  of  the  shock 
layer  and  A  from  the  solution. 

Initially  a  value  of  ~  0.1  was  chosen  based  on  the  order  of 

magnitude  study.  Once  this  ratio  was  fixed  we  had  twelve  equations  with 

twelve  unknowns  which,  in  theory  at  least,  are  soluble.  As  a  consequence, 

the  solution  failed  to  reach  the  hypothetical  body  located  at  a  known 

distance  from  the  shock.  In  order  to  surmount  this  difficulty,  the  shock 

standoff  distance  was  treated  as  an  additional  dependent  variable.  A 

comparison  of  values  of  is  shown  in  Figs.  11,  12  and  13.  It  was  found 

that  our  initial  guess  of  ■£—  ~  0.1  gave  too  large  a  standoff  distance.  On 

KB 

the  other  hand,  the  second  guess  *'■’  0.06  resulted  in  a  much  smaller 

value  of  A  .  The  actual  value  of  will  lie  between  these  two  limits 

provided  the  solution  converges. 

Since  a  solution  which  satisfies  the  boundary  conditions  for  the 
chemically  frozen  flow  on  the  stagnation  streamline  was  not  obtained,  no 
attempt  was  made  to  solve  the  reacting  case. 
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V.  CONCLUSIONS  AND  RECOMMENDATIONS 


Even  though  no  valid  solution  was  obtained  we  can  conclude  the 
following: 

1.  The  Newton-Raphson  iterative  method  is  applicable  for 
problems  of  this  kind  provided  the  initial  estimates  of  the  functions  and 
their  derivatives  are  fairly  good.  As  mentioned  in  the  text,  if  the  initial 
estimates  are  not  close  to  the  actual  values  no  convergence  will  be 
reached. 

2.  In  numerically  computing  the  elements  of  the  Jacobian,  one 
cannot  neglect  the  second  order  terms.  If  they  are  neglected  the  problem 
will  fail  to  converge. 

3.  The  value  of  the  shock  standoff  distance  cannot  be  specified  a 
priori.  It  has  to  be  obtained  as  a  part  of  the  solution. 

It  is  recommended  that  future  work  be  carried  on  in  an  attempt  to 
solve  the  whole  problem.  This  should  be  done  step  by  step  utilizing  the 
equations  in  a  nondimens ional  form.  First,  the  stagnation  streamline 
should  be  solved  for  a  frozen  case  and  then  following  the  same  method 
attempts  should  be  made  to  solve  the  reacting  case  with  pertinent  reaction 
equations.  The  solution  of  the  stagnation  streamline  will  constitute  an 
additional  boundary  from  which  the  iteration  for  the  whole  flow  field  can 
be  attempted. 

*! 
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Fig.  4  Values  Of  Functions  Of  Initial  Estimates  Vs.  r/Rg. 
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Fig.  5  Effect  Of  Initial  Slopes  On  Velocity  Ratio. 
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tance  0.06  Rg  And  0.10  Rg 


APPENDIX 


NONDIMENSIONAL  EQUATIONS  FOR  FROZEN  FLOW 
ON  THE  STAGNATION  LINE 


Conservation  of  Mass 


F.-rv£+2,V  +  ,y#  V-ff  =0 


Conservation  of  r-Momentum 


f2-T00Nr|/,u  +  T  dr'  J  2tT  a?;T+  dr1  dr'  J 

-ILL'^Ls^L  +  QL  +  4r.  J^_4u4.0 

r *\  dPle  5dQ+d6***r  dn  4u  J  u 


Conservation  of  Energy 

F>  ■  T- *  fe 'V  " 


j2i:  iii. 

wilY 

r'  dr' 

where 


REYNOLDS  NUMBER  NR=  R& 

PRANDTL  NUMBER  Np=  -fo  ^ 


koo 


dr'  W  dr'  p  dr' 


R0=  -gp  GAS  CONSTANT  FOR  THIS  GAS 
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* 

r'.J-  u'=-y- 

Rb  Uoo 

v'  =  JL. 
Uoo 

pi  = 

r  Pod 

T'  =  T/(  U®  ^ 

V*V 

II 

K8 

u/(£) 

\  cCpQ)/ 

and  the  9 -derivatives 

d\i'  Rell-rO 

are  assumed  to  have  the  forms: 

aV  _  -Rb 

d9  Rb  ■  Rs 

a^aa 

Rb~Rs 

ifll  Cpttt  /  J_  \  2Rb(I-'') 
dd 8  '  TTp  v  ^  "  /  Rb-Rs 


aV  .  RaO-rM 
a©2  "  ^(Rb-Rs) 
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